The purpose of this Python notebook is twofold.
The first is to serve as a quick start guide where you should be able to get started with the package by simply looking at the examples here and copying them to your liking.
The second is as a unit testing replacement. It is hard to unit test stochastic algorithms as the output may not (and should not) be the same thing every time. Therefore, instead, if all the examples included below work well, then you can assume that the package installed correctly and is working fine.
Before, getting started, we start by doing some basic plotting configuration and importing the numpy library. Advanced users can modify this to their liking.
In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
import matplotlib as mpl
#%config InlineBackend.figure_f.ormats=['svg']
color_list = ['r', 'k', 'b','g','y','m','c']
mpl.rc('axes', prop_cycle=(mpl.cycler('color', color_list) ))
mpl.rc('xtick', labelsize=12)
mpl.rc('ytick', labelsize=12)
import numpy as np
Bioscrape models consist of a set of species and a set of reactions (delays will be discussed later). These models can be simulated either stochastically via SSA or deterministically as an ODE. Each reaction is of the form
$\{INPUTS\} \xrightarrow[]{\rho(.)} \{OUTPUTS\}$
Here, INPUTS represent a multiset of input species and OUTPUTS represents a multiset of output species. The function $\rho(.)$ is either a deterministic rate function or a stochastic propensity. Propensities are identified by their name and require parameter dictionary with the appropriate parameters. The following functions are supported:
More details on all these propensity types can be found in the wiki documentation
First, the following model of transcription and translation will be created programatically. There are three chemical species: $G$ is a gene, $T$ is a transcript, $X$ is a protein.
The first reaction uses a proportional positive hill function as its rate function to represent induction. The second reaction uses a positive hill function function to represent ribosome saturation. The third and fourth reactions reaction represent degredation via dilution. No delays will be included in this model. This model is constructed below and simulated both stochastically and deterministically.
In [2]:
from bioscrape.simulator import py_simulate_model
from bioscrape.types import Model
#Create a list of species names (strings)
species = ["G", "T", "X", "I"]
#create a list of parameters in the form (param_name[string], param_val[number])
params = [("ktx", 1.5), ("ktl", 10.0), ("KI", 10), ("n", 2.0), ("KR", 20), ("delta", .1)]
#create reaction tuples in the form:
#(Inputs[string list], Outputs[string list], propensity_type[string], propensity_dict {propensity_param:model_param})
rxn1 = (["G"], ["G", "T"], "proportionalhillpositive", {"d":"G", "s1":"I", "k":"ktx", "K":"KI", "n":"n"})
rxn2 = (["T"], ["T", "X"], "hillpositive", {"s1":"T", "k":"ktl", "K":"KR", "n":1}) #Notice that parameters can also take numerical values instead of being named directly
rxn3 = (["T"], [], "massaction", {"k":"delta"})
rxn4 = (["X"], [], "massaction", {"k":"delta"})
#Create a list of all reactions
rxns = [rxn1, rxn2, rxn3, rxn4]
#create an initial condition dictionary species not included in the dictionary will default to 0
x0 = {"G":1, "I":10}
#Instaniate the Model object
M = Model(species = species, parameters = params, reactions = rxns, initial_condition_dict = x0)
#Simulate the Model deterministically
timepoints = np.arange(0, 150, .1)
results_det = py_simulate_model(timepoints, Model = M) #Returns a Pandas DataFrame
#Simulate the Model Stochastically
results_stoch = py_simulate_model(timepoints, Model = M, stochastic = True)
#Plot the results
plt.figure(figsize = (12, 4))
plt.subplot(121)
plt.title("Transcript T")
plt.plot(timepoints, results_det["T"], label = "deterministic")
plt.plot(timepoints, results_stoch["T"], label = "stochastic")
plt.legend()
plt.xlabel("Time")
plt.subplot(122)
plt.title("Protein X")
plt.plot(timepoints, results_det["X"], label = "deterministic")
plt.plot(timepoints, results_stoch["X"], label = "stochastic")
plt.legend()
plt.xlabel("Time")
Out[2]:
In stochastic simulations, bioscrape also supports delay. In a delay reaction, delay inputs/outputs are consumed/produced after some amount of delay. Reactions may have a mix of delay and non-delay inputs and outputs. Bioscrape innately supports a number of delay-types:
In the following example model, the following delays are added to the transcription and translation reactions described above and then simulated stochastically. Note that delays and delay inputs/outputs will be ignored if a model with delays is simulated deterministically.
In [3]:
from bioscrape.simulator import py_simulate_model
from bioscrape.types import Model
#create reaction tuples with delays require additional elements. They are of the form:
#(Inputs[string list], Outputs[string list], propensity_type[string], propensity_dict {propensity_param:model_param},
# delay_type[string], DelayInputs [string list], DelayOutputs [string list], delay_param_dictionary {delay_param:model_param}).
rxn1d = (["G"], ["G"], "proportionalhillpositive", {"d":"G", "s1":"I", "k":"ktx", "K":"KI", "n":"n"},
"gaussian", [], ["T"], {"mean":10.0, "std":1.0})
rxn2d = (["T"], ["T"], "hillpositive", {"s1":"T", "k":"ktl", "K":"KR", "n":1},
"gamma", [], ["X"], {"k":10.0, "theta":3.0})
#Reactions 3 and 4 remain unchanged
rxns_delay = [rxn1d, rxn2d, rxn3, rxn4]
#Instaniate the Model object, species, params, and x0 remain unchanged from the previous example
M_delay = Model(species = species, parameters = params, reactions = rxns_delay, initial_condition_dict = x0)
#Simulate the Model with delay
results_delay = py_simulate_model(timepoints, Model = M_delay, stochastic = True, delay = True)
#Plot the results
plt.figure(figsize = (12, 4))
plt.subplot(121)
plt.title("Transcript T")
plt.plot(timepoints, results_det["T"], label = "deterministic (no delay)")
plt.plot(timepoints, results_stoch["T"], label = "stochastic (no delay)")
plt.plot(timepoints, results_delay["T"], label = "stochastic (with delay)")
plt.legend()
plt.xlabel("Time")
plt.subplot(122)
plt.title("Protein X")
plt.plot(timepoints, results_det["X"], label = "deterministic (no delay)")
plt.plot(timepoints, results_stoch["X"], label = "stochastic (no delay)")
plt.plot(timepoints, results_delay["X"], label = "stochastic (with delay)")
plt.legend()
plt.xlabel("Time")
Out[3]:
In deterministic and stochastic simulations, bioscrape also supports rules which can be used to set species or parameter values during the simulation. Rules are updated every simulation timepoint - and therefore the model may be sensitive to how the timepoint spacing.
The following example two rules will be added to the above model (without delay).
Rules can also be used for quasi-steady-state or quasi-equilibrium approximations, computing parameters on during the simulation, and much more!
There are two main types of rules:
In [4]:
#Create new parameters for rule 1. Model is now being modified
M.create_parameter("I0", 10) #Inducer concentration
M.create_parameter("T", 25) #Initial time inducer is added
#Create rule 1:
#NOTE Rules can also be passed into the Model constructor as a list of tuples [("rule_type", {"equation":"eq string"})]
M.create_rule("assignment", {"equation":"I = _I0*Heaviside(t-_T)"}) #"_" must be placed before param names, but not species.
#Rule 2 will use constants in equations instead of new parameters.
M.create_rule("assignment", {"equation":"S = 50*X/(1+.2*X)"}) #Species S is added automatically
#reset the initial concentration of the inducer to 0
M.set_species({"I":0})
#Simulate the Model deterministically
timepoints = np.arange(0, 150, 1.0)
results_det = py_simulate_model(timepoints, Model = M) #Returns a Pandas DataFrame
#Simulate the Model Stochastically
results_stoch = py_simulate_model(timepoints, Model = M, stochastic = True)
#Plot the results
plt.figure(figsize = (12, 8))
plt.subplot(223)
plt.title("Transcript T")
plt.plot(timepoints, results_det["T"], label = "deterministic")
plt.plot(timepoints, results_stoch["T"], label = "stochastic")
plt.legend()
plt.subplot(224)
plt.title("Protein X")
plt.plot(timepoints, results_det["X"], label = "deterministic")
plt.plot(timepoints, results_stoch["X"], label = "stochastic")
plt.legend()
plt.subplot(221)
plt.title("Inducer I")
plt.plot(timepoints, results_det["I"], label = "deterministic")
plt.plot(timepoints, results_stoch["I"], label = "stochastic")
plt.legend()
plt.xlabel("Time")
plt.subplot(222)
plt.title("Signal S")
plt.plot(timepoints, results_det["S"], label = "deterministic")
plt.plot(timepoints, results_stoch["S"], label = "stochastic")
plt.legend()
plt.xlabel("Time")
Out[4]:
In [5]:
M.write_bioscrape_xml('models/txtl_model.xml')
f = open('models/txtl_model.xml')
print("Bioscrape Model XML:\n", f.read())
M_loaded = Model('models/txtl_model.xml')
#Change the induction time
#NOTE That changing a model loaded from xml will not change the underlying XML.
M_loaded.set_parameter("T", 50)
#Simulate the Model deterministically
timepoints = np.arange(0, 150, 1.0)
results_det = py_simulate_model(timepoints, Model = M_loaded) #Returns a Pandas DataFrame
#Simulate the Model Stochastically
results_stoch = py_simulate_model(timepoints, Model = M_loaded, stochastic = True)
#Plot the results
plt.figure(figsize = (12, 8))
plt.subplot(223)
plt.title("Transcript T")
plt.plot(timepoints, results_det["T"], label = "deterministic")
plt.plot(timepoints, results_stoch["T"], label = "stochastic")
plt.legend()
plt.subplot(224)
plt.title("Protein X")
plt.plot(timepoints, results_det["X"], label = "deterministic")
plt.plot(timepoints, results_stoch["X"], label = "stochastic")
plt.legend()
plt.subplot(221)
plt.title("Inducer I")
plt.plot(timepoints, results_det["I"], label = "deterministic")
plt.plot(timepoints, results_stoch["I"], label = "stochastic")
plt.legend()
plt.xlabel("Time")
plt.subplot(222)
plt.title("Signal S")
plt.plot(timepoints, results_det["S"], label = "deterministic")
plt.plot(timepoints, results_stoch["S"], label = "stochastic")
plt.legend()
plt.xlabel("Time")
Out[5]:
The next cell imports a model from an SBML file and then simulates it using a deterministic simulation. There are limitations to SBML compatibility.
Below, we first plot out the simulation results for an SBML model where a species X0 goes to a final species X1 through an enymatic process.
In [6]:
from bioscrape.sbmlutil import import_sbml
M_sbml = import_sbml('models/sbml_test.xml')
timepoints = np.linspace(0,100,1000)
result = py_simulate_model(timepoints, Model = M_sbml)
plt.figure()
for s in M_sbml.get_species_list():
plt.plot(timepoints, result[s], label = s)
plt.legend()
Out[6]:
We plot out the repressilator model found here. This model generates oscillations as expected. Highlighting the utility of this package, we then with a single line of code switch to a stochastic simulation and note that the amplitudes of each burst become noisy.
In [7]:
# Repressilator deterministic example
import bioscrape
M_represillator = bioscrape.sbmlutil.import_sbml('models/repressilator_sbml.xml')
#Simulate Deterministically and Stochastically
timepoints = np.linspace(0,700,10000)
result_det = py_simulate_model(timepoints, Model = M_represillator)
result_stoch = py_simulate_model(timepoints, Model = M_represillator, stochastic = True)
#Plot Results
plt.figure(figsize = (12, 8))
for i in range(len(M_represillator.get_species_list())):
s = M_represillator.get_species_list()[i]
plt.plot(timepoints, result_det[s], color = color_list[i], label = "Deterministic "+s)
plt.plot(timepoints, result_stoch[s], ":", color = color_list[i], label = "Stochastic "+s)
plt.title('Repressilator Model')
plt.xlabel('Time')
plt.ylabel('Amount')
plt.legend()
Out[7]:
In [ ]: